home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / xvecmath.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  7KB  |  304 lines
  1. /*
  2.  *    Math definitions
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)xvecmath.cc    2.1    8/18/89
  23.  */
  24.  
  25. #include "rw/<T>Vec.h"
  26. #define TYPE <T>_TYPE
  27. #include "vecdefs.h"
  28.  
  29. // The type returned by abs depends on the type of vector
  30. #if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
  31. #define ABS_TYPE <P>
  32. #define ABS_VEC <P>Vec
  33. #else
  34. #define ABS_TYPE <T>
  35. #define ABS_VEC <T>Vec
  36. #endif
  37.  
  38. #ifdef APPLY_MATH
  39. #if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
  40. <T>Vec
  41. <T>Vec::apply(CmathFunTy f)
  42. // Apply function returning <T> to each element of vector slice
  43. {
  44.   register i = length();
  45.   <T>Vec temp(i);
  46.   register <T>* sp = data();
  47.   register <T>* dp = temp.data();
  48.   register j = stride();
  49.   while (i--) { *dp++ = f(*sp);  sp += j; }
  50.   return temp;
  51. }
  52.  
  53. <P>Vec
  54. <T>Vec::apply2(CmathFunTy2 f)
  55. // Apply function returning <P> to each element of vector slice
  56. {
  57.   register i = length();
  58.   <P>Vec temp(i);
  59.   register <T>* sp = data();
  60.   register <P>* dp = temp.data();
  61.   register j = stride();
  62.   while (i--) { *dp++ = f(*sp);  sp += j; }
  63.   return temp;
  64. }
  65.  
  66. // Version 2.0 has trouble finding the proper overloaded real & imag.
  67. // Put in an explicit version:
  68. <P>Vec
  69. real(const <T>Vec& v)
  70. {
  71.   register i = v.length();
  72.   <P>Vec temp(i);
  73.   register <T>* sp = v.data();
  74.   register <P>* dp = temp.data();
  75.   register j = v.stride();
  76.   while(i--) {*dp++ = real(*sp); sp += j;}
  77.   return temp;
  78. }
  79.  
  80. <P>Vec
  81. imag(const <T>Vec& v)
  82. {
  83.   register i = v.length();
  84.   <P>Vec temp(i);
  85.   register <T>* sp = v.data();
  86.   register <P>* dp = temp.data();
  87.   register j = v.stride();
  88.   while(i--) {*dp++ = imag(*sp); sp += j;}
  89.   return temp;
  90. }
  91.  
  92. #else  // Not Complex
  93. <T>Vec
  94. <T>Vec::apply(mathFunTy f)
  95. // Apply function returning <T> to each element of vector slice
  96. {
  97.   register i = length();
  98.   <T>Vec temp(i);
  99.   register <T>* sp = data();
  100.   register <T>* dp = temp.data();
  101.   register j = stride();
  102.   while (i--) { *dp++ = f(*sp);  sp += j; }
  103.   return temp;
  104. }
  105. #endif
  106. #endif
  107.  
  108. ABS_VEC
  109. abs(const <T>Vec& s)
  110. // Absolute value of <T>Vec slice.
  111. {
  112.   register i = s.length();
  113.   ABS_VEC temp(i);
  114.   register <T>* sp = s.data();
  115.   register ABS_TYPE* dp = temp.data();
  116.   register j = s.stride();
  117.   while (i--) { *dp++ = ABSOLUTE(*sp);  sp += j; }
  118.   return temp;
  119. }
  120.  
  121. #if HAS_ATAN2
  122. <T>Vec
  123. atan2(const <T>Vec& u,const <T>Vec& v)
  124.      // Arctangent of u/v.
  125. {
  126.   register i = v.length();
  127.   u.lengthCheck(i);
  128.   <T>Vec temp(i);
  129.   register <T>* up = u.data();
  130.   register <T>* vp = v.data();
  131.   register <T>* dp = temp.data();
  132.   register uj = u.stride();
  133.   register vj = v.stride();
  134.   while (i--) { *dp++ = atan2(*up,*vp);  up += uj;  vp += vj; }
  135.   return temp;
  136. }
  137. #endif
  138.  
  139. <T>Vec
  140. cumsum(const <T>Vec& s)
  141. // Cumulative sum of <T>Vec slice.
  142. // Note: V == delta(cumsum(V)) == cumsum(delta(V))
  143. {
  144.   unsigned i = s.length();
  145.   <T>Vec temp(i);
  146.   register <T>* sp = s.data();
  147.   register <T>* dp = temp.data();
  148.   REGISTER <T> c = 0;
  149.   register j = s.stride();
  150.   while (i--) { c += *sp;  *dp++ = c;  sp += j; }
  151.   return temp;
  152. }
  153.  
  154. <T>Vec
  155. delta(const <T>Vec& s)
  156.      // Element differences of <T>Vec slice.
  157.      // Note: V == delta(cumsum(V)) == cumsum(delta(V))
  158. {
  159.   unsigned i = s.length();
  160.   <T>Vec temp(i);
  161.   register <T>* sp = s.data();
  162.   register <T>* dp = temp.data();
  163.   REGISTER <T> c = 0;
  164.   register j = s.stride();
  165.   while (i--) { *dp++ = *sp - c;  c = *sp;  sp += j; }
  166.   return temp;
  167. }
  168.  
  169. <T>
  170. dot(const <T>Vec& u,const <T>Vec& v)
  171.      // Vector dot product.
  172. {
  173.   register i = v.length();
  174.   u.lengthCheck(i);
  175.   REGISTER <T> t =0;
  176.   register <T>* up = u.data();
  177.   register <T>* vp = v.data();
  178.   register uj = u.stride();
  179.   register vj = v.stride();
  180.   while (i--) { t += *up * *vp;  up += uj;  vp += vj; }
  181.   return t;
  182. }
  183.  
  184. #if HAS_MINMAX
  185. int
  186. max(const <T>Vec& s)
  187. // Index of maximum element.
  188. {
  189.   if (s.length() == 0) s.emptyErr("max");
  190.   register <T>* sp = s.data();
  191.   REGISTER <T> t = *sp;
  192.   register j = s.stride();
  193.   register x = 0;
  194.   for (register i =0; i<s.length(); i++) {
  195.     if (*sp > t) { t = *sp;  x = i; }
  196.     sp += j;
  197.   }
  198.   return x;
  199. }
  200.  
  201. int
  202. min(const <T>Vec& s)
  203.      // Index of minimum element.
  204. {
  205.   if (s.length() == 0) s.emptyErr("min");
  206.   register <T>* sp = s.data();
  207.   REGISTER <T> t = *sp;
  208.   register j = s.stride();
  209.   register x = 0;
  210.   for (register i =0; i<s.length(); i++) {
  211.     if (*sp < t) { t = *sp;  x = i; }
  212.     sp += j;
  213.   }
  214.   return x;
  215. }
  216. #endif
  217.  
  218. <T>
  219. prod(const <T>Vec& s)
  220.      // Product of <T>Vec elements.
  221. {
  222.   register i = s.length();
  223.   register <T>* sp = s.data();
  224. #if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
  225.   REGISTER <T> t = <T>(1,1);
  226. #else
  227.   REGISTER <T> t = 1;
  228. #endif
  229.   register j = s.stride();
  230.   while (i--) { t *= *sp;  sp += j; }
  231.   return t;
  232. }
  233.  
  234. #if HAS_POW
  235. <T>Vec
  236. pow(const <T>Vec& u,const <T>Vec& v)
  237. // u to the v power.
  238. {
  239.   register i = v.length();
  240.   u.lengthCheck(i);
  241.   <T>Vec temp(i);
  242.   register <T>* up = u.data();
  243.   register <T>* vp = v.data();
  244.   register <T>* dp = temp.data();
  245.   register uj = u.stride();
  246.   register vj = v.stride();
  247.   while (i--) { *dp++ = pow(*up,*vp);  up += uj;  vp += vj; }
  248.   return temp;
  249. }
  250. #endif
  251.  
  252. <T>Vec
  253. reverse(const <T>Vec& s)
  254.      // Reverse vector elements.
  255. {
  256.   register i = s.length();
  257.   <T>Vec temp(i);
  258.   register <T>* sp = s.data();
  259.   register <T>* dp = &temp(i-1);
  260.   register j = s.stride();
  261.   while (i--) { *(dp--) = *sp;  sp += j; }
  262.   return temp;
  263. }
  264.  
  265. <T>
  266. sum(const <T>Vec& s)
  267.      // Sum of a <T>Vec
  268. {
  269.   register int i = s.length();
  270.   register <T>* sp = s.data();
  271.   REGISTER <T> t = 0;
  272.   register j = s.stride();
  273.   while (i--) { t += *sp;  sp += j; }
  274.   return t;
  275. }
  276.  
  277. #if HAS_VARIANCE
  278. <P>
  279. variance(const <T>Vec& s)
  280. // Variance of a <T>Vec
  281. {
  282.   register int i = s.length();
  283.   if(i<2){
  284.     RWnote("::variance()", "Less than two points.");
  285.     RWerror(DEFAULT);
  286.   }
  287.   register <T>* sp = s.data();
  288.   register <P> sum2 = 0;
  289.   REGISTER <T> avg = mean(s);
  290.   register j = s.stride();
  291.   while (i--) {
  292.     REGISTER <T> temp = *sp - avg;
  293. #if TYPE==DComplex_TYPE || TYPE==FComplex_TYPE
  294.     sum2 += norm(temp);
  295. #else
  296.     sum2 += temp * temp;
  297. #endif
  298.     sp += j;
  299.   }
  300.   // Use the biased estimate (easier to work with):
  301.   return sum2/s.length();
  302. }
  303. #endif
  304.